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ABSTRACT 

We use the high-resolution Aquarius simulations of the formation of Milky Way- 
sized haloes in the ACDM cosmology to study the effects of dark matter substruc- 
tures on gravitational lensing. Each halo is resolved with ~ 10^ particles (at a mass 
resolution irip ~ 10'^ to 10^/i~^Mo) within its virial radius. Subhaloes with masses 
TOsub ^ lO^/i~^M0 are well resolved, an improvement of at least two orders of mag- 
nitude over previous lensing studies. We incorporate a baryonic component modelled 
as a Hernquist profile and account for the response of the dark matter via adiabatic 
contraction. We focus on the "anomalous" flux ratio problem, in particular on the vi- 
olation of the cusp-caustic relation due to substructures. We find that subhaloes with 
masses less than ~ IO^H'^Mq play an important role in causing flux anomalies; such 
low mass subhaloes have been unresolved in previous studies. There is large scatter 
in the predicted flux ratios between different haloes and between different projections 
of the same halo. In some cases, the frequency of predicted anomalous flux ratios is 
comparable to that observed for the radio lenses, although in most cases it is not. 
The probability for the simulations to reproduce the observed violations of the cusp 
lenses is w 10~^. We therefore conclude that the amount of substructure in the cen- 
tral regions of the Aquarius haloes is insufficient to explain the observed frequency of 
violations of the cusp-caustic relation. These conclusions are based purely on our dark 
matter simulations which ignore the effect of baryons on subhalo survivability. 

Key words: Gravitational lensing - dark matter - galaxies: ellipticals - galaxies: 
formation 



1 INTRODUCTION 

Currently there are ~ 200 known galaxy-scale lenses, di- 
vided roughly equally in number into lensed active galac- 
tic n uclep and lensed background galaxies (|Bolton et al.l 
I2OO8I ). These galaxy-scale lenses allow diverse applications 
(see the review "Strong Gravitational Lensing" by Kochanek 
in[Schnoidcr ct al. 2006) such as a determination of the Hub- 
ble constant, a characterisation of galaxy evolution, and 
measurements of the mass distribution in galaxies. The last 
application will likely be the most important one in the next 
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decade, since there are few other probes at intermediate red- 
shifts (2 ~ 0.5 - 1). 

It was noticed quite early on that the flux ratios of 
the multi-lensed images are more difficult to reproduce with 
simple parametr ic mass models than the image positions 
llKochanekHigQll l. This has been termed the "anomalous flux 
ratio" problem. Image positions and magnifications (flux 
ratios) are determined by the flrst-order and second-order 
derivatives of the lensing potential, respectively. Therefore 
flux ratios, as a high-order derivative, are expected to be 
more sensitive to small changes in the lensing potential than 
image positions. 

In this regard, gravitational lenses with two or three 
close images deserve special attention because, in these 
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cases, the sources must be close to either a fold or a cusp 
of the caustic. It is well known for any smooth lensing 
potential that the close images follow asymptotic flux ra- 
tio relations: for a close pair, their flux ratio approaches 
unity when their separation goes to zero, while for a close 
triple, the ratio of the flux of the middle image to the 
sum of the fluxes o f the two outer images asymptoti- 
cally goes to unity (iMaol Il992!: 'Schneider fc WeissI (l992l : 
iKeeton et al.ll2003l : ICongdon et al. 2008). However, the ob- 
served lensing systems often violate these asymptotic re- 
lations. This wa s taken to be evid ence for substructure in 
lensing galaxies (iMao fc Schneider|[l99 8: Motcalf fc Madau 
_2001; M etcalf fc Z hao 2002; Dal ai fc Kochanekll2002i : IChibal 
2002. ; .Kochanek fc Dalai 2004 ) on the physical scale of 
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the separation between close images (typically of the 
order of ~ 1 kpc). Spectroscopic observations go be- 
yond simple broad-band flux ratios and pro vide a promis- 
ing way to probe sub s tructure in lenses (|Metcalf et al.l 
l2004l : IChiba et al.l l2005l : ISugai et al.l bOOTh . Other sugges- 
tive evidence for sub structures comes from astrometry 
(see IChen et al.' '20071 for a general discussion), such as 
bent jets (.Metcalf , 200 j). and deta il ed image structures for 
B2016-I-112 JSchneider et al.l liooil: iKoopmans et all |2002|: 



More et al.ll2009D and B0128+437 (|Biggs et al.ll2004l : Izhan3 
2008h . Substructures may also have detectable effects on 



the t ime delays in gravitational lenses (|Keeton fc MoustakasI 

|200^ 

lEvans fc Wit j (gooi) argued that some of these lensing 
"anomahes" may be accommodated by changes in the po- 
tentials of the main lensing galaxies in parametric models. 
However, signi ficant changes are needed in order to explain 
the a nomalies (|Kochanek fc Dalalll2004l : ICongdon fc KeetonI 
l2005l ). The angular structures of the lenses whenever the 
measurements are available suggest ellipsoidal central poten- 
tials, where the high amplitude, higher order multipoles that 
are required to explain the flux ratio anom alies arc not seen 
IIKochanek fc Dalall 12004 lYoo et al. I l2005l , 12006 ). There are 
strong hints that substructures may indeed be real in lens- 
ing galaxies. First, observationally, saddle (negative-parity) 
images are often fainter than the predictions of simple 
smooth models. This is expected from lensing by substruc- 
ture (such as stars, or subhalo es; ISchechter fc WambsganssI 
l2002l :l Kochanek fc Dalaill2004l ) , but impossible to explain by 
propagation effects, such as galactic scintillation and scatter 
broadening, as earlier postulated (jKoopma ns ct al. 2003). 
This arguably constitutes the most convincing evidence for 
substructure lensing. Second, in many gravitational lenses, 
the substructure is directly seen as luminous satellites. For 
example , nearly half of the CLASS lenses (|Browne et al.l 
[2003 ; Myers etal.ll2003l : I Jackson et al.ll2009l ) show luminous 
satellite galaxies within a few kpc of the primary lensing 
galaxie^. Inclusion of satellites in the modelling dramati- 
cally improves the fit to the image positions. In the case of 
B2045+265, the inclusion of a com panion galaxy he lps to 
explain the flux ratio anomaly (.McKean et al.l [20071 ). The 
additional dark subhaloes within the main lensing galaxies, 
as well as the intergalactic perturbers along the line-of-sight 



This fract i on is a factor of ~ 2 higher than that claimed in 
I Bryan et al.l ||2008| ) as reveale d by a more careful analysis of HST 
images of the CLASS lenses ll Jackson et al.ir2009l) . 



[Miranda fc Macciol 20071 ) mav also help to explain the ob- 
served lensing anomalies. 

Much of the interest in (milli-) lensing flux anomalies 
arises because they may be caused by the elusive sub- 
structure generically predicted by the hierarchical struc- 
ture formation in the cold dark matter (CDM) cosmology 
(e.g.[ Kauffmann ct al. 1993; Klvpin ot al. 1999; Moor e et al 
1999ljGhigna et al...2000. ; ,Gao et al...2004 a,b; .Diemand et al 
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20071 ) ■ In this model, large structures form via merging 
and accretion of smaller structures. The cores of these 
small structures often survive tidal destruction and man- 
ifest themselves as subhaloes (substructure). Recent high- 
resolution simulations predict many thousands of subhaloes 
(down to rriau b ~ 10^ Mq, or to circular velocity of V c ~ 
4kms-^; e.g. iMadau et all [20081 : ISpringel et al.ll2008l ), at 
least two orders of magnitude more than the number of ob- 
served satellite galaxies in the Milky Way, even after ac- 
counting for the newly discovere d faint satellite galax ies 
from the Sloan Digital Sky Survey ([Belokurov et al.ll2007l ). A 
possible solution is that star formation may be strongly sup- 



Efstathioulll992l:lKauffmann et al.lll99a 


: iThoul fc Weinberd 


199& .Bullock et al. 2000l: Gnedin|[2000l: 


Benson et al.[ 20021), 



and thus they remain dark and difficult to detect through 
light-based methods. If this is the case, then gravitational 
lensing can potentially probe this population since it de- 
pends only on the mass but not on whether the lenses are 
luminous or dark. 

Numerical simulations indicate that subhaloes typi- 
cally account for 5- 10 per cent of the total mass in a 
galaxy-type halo (e .g. iKlvpin et al.lll999l:lMoore et al.lll99£ 



[Ghigna et a"l][2000l ) . The study bv [Dalal fc Kochanekl ([2002 ) 

requires fsuh ~ 0.6% to 7% (with a median of 2%) of the 
mass to be in substructures (90% confidence limit) in or- 
der to explain the observed fiux anomaly problem. At first 
sight, the fraction of substructure from simulations seems to 
be more than sufficient to explain the flux anomaly. Upon 
closer examination, however, a problem emerges: lensing 
probes the central few kpc around the line-of-sight through 
the galaxy, while most substructures are in the outer re- 
gions of its dark matter halo, since those that come close 
to the centre are tidally destroyed. Thus it remains un- 
clear whether the predicted substructure in the inner regions 
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). In contrast, on cluster scales, the 



amount of predicted substructure seems to be consisten t 
with weak and strong lensing data ([Nataraian et al]|2007l ). 

Previous lensing studies simulated galaxy-sized haloes 
with ~ 10® particles so that subhaloes were resolved down to 
~ lO'^ to 10^ h~ ^ Mq. State-of-the-art simulations can now 
resolve haloes with two or even three orders of magnitude 
more particles, thus reaching substantially lower mass sub- 
haloes. In this work, we revisit the issue of substructure lens- 
ing using the Aquarius simulations of six galaxy-sized haloes. 
These coUisionless A'^-body simulations were performed by 
the Virgo Consortium in a concordance ACDM universe. 
The subhaloes in ea ch halo are resolved d own to masses of 
m-sub ~ 10^ h~^MQ ([Springel et al.ll2008l) , at least two or- 
ders of magnitude better than that in previous substructure 
lensing studies. 
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Our paper is organised as follows. In Section 2, we de- 
scribe the realisation and the properties of the simulated 
lensing galaxies. Our methods and techniques for the lens- 
ing simulations together with our test results are presented 
in Section 3. In Section 4, we apply our lensing simulation 
to the six simulated galaxy haloes from the Aquarius simu- 
lation to derive their lensing properties, including the cusp 
relations, and we compare the numerical results with ob- 
servations. A summary of the paper and a discussion are 
given in Section 5. The cosmology we adopt for the lensing 
simulatio n is the same as tha t used for the Aquarius sim- 
ulations ijSpringel et al ] |2005l ). with a matter density ^^ni 
= 0.25, cosmological constant f^A ~ 0.75, Hubble constant 
h = _ffo/(100 kms~^ Mpc~^) — 0.73 and linear fluctuation 
amplitude erg — 0.9. 
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2 FROM DARK MATTER HALOES TO 
EARLY- TYPE LENSING GALAXIES 

In this section, we summarise the properties of dark matter 
haloes from the Aquarius simulations relevant to our study, 
in particular the subh alo properties. Readers are referred to 
ISprineel et al] (|2008l ) for more details. We will show that 
dark matter alone is, as expected, insufficient to cause mul- 
tiple image splittings, and therefore we must incorporate a 
stellar component; we detail such a procedure in ^2.2\ 



2.1 The Aquarius simulations 

The Aquarius project (jSpringel et al.|[2008l ) is a suite of sim- 
ulations of six galaxy-sized dark matter haloes with five lev- 
els of numerical resolution. The haloes were selected from a 
100h~^ Mpc simulation box within the concordance cosmol- 
ogy (for parameters see above). The simulations were run 
with GADGET-3, an improved version o f the GADGET-2 
code ijSpringel et al.ll200ll : ISpringelll2005h . The highest res- 
olution level (level 1) was achieved for only one halo {"Aq- 
A-1") with ~ 1.5 billion halo particles. Level-2 simulations 
were performed for a sample of six dark matter haloes, with 
about 200 million particles per halo. The softening length 
is ~0.05/i~^ kpc, and the mass resolution ranges from 10^ 
to lO''/i~^M0. All haloes are Milky- Way type systems in 
terms of their mass and rotation curve. We will use the six 
level-2 haloes {Aq-A-2, Aq-B-2, Aq-C-2, Aq-D-2, Aq-E-2 and 
Aq-F-2) at redshift zero for our analysis of substructure lens- 
ing. As we will show later on, the scatter in lensing proper- 
ties among different haloes (and for different projections) is 
large, and so it is important to examine more than one halo 
for statistical purposes. 

The basic properties of the six haloes at z = are listed 
in Table [U In particular, all the density profiles are rea- 
sonably fit by Navarro, Frenk and White (NEW) profiles 



Figure 1. Density profiles (solid curves), multiplied by r^, for 
the halo Aq-D-2 at redshifts 2 = 0, 0.5 and 0.99. All haloes are 
reasonably well fitted by the NFW profile (dotted curves, see 
Eq. [T}, which follows p(r) oc on small scales and p{r) oc 
on large scales. The vertical dashed lines indicate the softening 
length and r2oo- 

(|Navarro et al.ll 19961 . Il997l fl: 

^ M200 

4Tvr{r + r 200/ cyfic)' 
M(-:r) - M20of{rc/r2oo) (1) 
^ ^ f(c) 
/(c) =ln(l + c)-c/(l + c), 

where r2oo is the radius within which the mean dark halo 
mass density is 200 times the critical density, M200 is the 
mass enclosed within 7-200 , and c = r2oo / ''s is the concentra- 
tion parameter with Ts being the scale radius. 

We artificially put all these haloes (snapshot z — 0.0) 
at redshift z = 0. 6 (corresponding ro ughly to the most likely 
lens redshift, e.g. [Turner et al]|l984l ). keeping their physical 
sizes unchanged. However, we also take a snapshot of the 
halo Aq-A-2 at redshift z — 0.6 as a lens, and compare its 
lensing properties with those artificially shifted to z = 0.6. 
As will be shown in §4, the scatter among the six haloes is 
much larger than the differences between haloes at redshifts 
z = and z — 0.6, and so adopting the z = haloes will not 
significantly change the properties of substructure lensing. 
This is also seen in the evolution of density profiles of these 
haloes. Eig.[T] shows the density profiles for the halo Aq-D- 
2 at redshifts 0, 0.50 and 0.99. The changes in the profiles 
since redshift 1 are relatively small since the Aquarius haloes 
form earlier than that. 

As we are primarily interested in the substructure lens- 
ing, an important step is the identification of the subhaloes. 

^ An even better fit is found using the lEinastol lll966l ') profile 
llNavarro et al.ir2008t) , but here we adopt the simpler NFW profile 
which we use later to take into account the adiabatic contraction 
of dark matter haloes. 



4 Xu et al. 



Table 1. Dark matter halo properties in the Aquarius simulations: 



Halo Name 


'•200 

kpc) 


Mtot 
(IQiO/i-iM©) 


c 


Mass Resolution 
(h-iMo) 


^^200 


A'sub 


/sub 

(per cent) 


Aq-A-2 


179.5 


132.8 


16.2 


1.0 X 10"* 


1.3 X 10* 


2.1 X 10** 


7.14 


Aq-B-2 


137.1 


59.5 


9.7 


4.7 X 10^ 


1.3 X 10* 


2.5 X 10* 


6.98 


Aq-C-2 


177.3 


127.7 


15.2 


1.0 X 10"* 


1.2 X 10* 


1.7 X 10* 


4.12 


Aq-D-2 


177.3 


128.5 


9.4 


1.0 X 10* 


1.3 X 10* 


2.2 X 10* 


6.56 


Aq-E-2 


155.0 


85.7 


8.3 


7.0 X 10^ 


1.2 X 10* 


2.3 X 10* 


7.28 


Aq-F-2 


153.0 


80.5 


9.8 


4.9 X 10^ 


1.6 X 10* 


2.6 X 10* 


11.20 


Aq-A-2 {Z = 0.6) 


134.4 


92.2 


10.4 


1.0 X 10"* 


9.3 X 10^ 


1.7 X 10* 


6.50 



Note: Col (1): halo name, Cols (2)-(4): r200i c and M200 are defined in eq. JTJl for the main halo. Mtot = jV/200 + ^sub, where A/tot and 
-'^sub E're the total masses of all dark matter and of all the subhaloes within r200. Col (5): Mass resolution {h~^ Mq). Col (6): Af200 is 
the total number of particles within r2oo- Col (7): A^aub is the number of subhaloes within r2oo. Col (8): fsnb is the mass fraction of 
subhaloes within r2oa, defined by A/gub/Mtot- 



We use the SUBFIND routine (|Springelll2005l ) to identify sub- 
haloes exceeding 20 particles, which corresponds to a min- 
imum subhalo mass of ~ 10^ h-'^MQ. The number of sub- 
haloes in each halo ranges from about 1.7 x 10'' to 2.6 x 10* 
within r2ao, with 4.1-11.2 per cent of the total halo mass 
locked up in bound subhaloes (see Col (8): faub in Table [T|. 

The subhalo mass f unction follows a power-law: 
dA'^(msub)/dmsub oc m^^j^® (|SpringeI et al.|[200^ . The aver- 
age mass of subhaloes (within r2oo) is ~ 10® to 10''^h~^MQ 
and their average half-mass radius is < 0.2h~^ kpc, with 
large scatter. The most massive subhalo has a mass of ~ 10^ 
to IO^^/i-^Mq and a half- mass radius ~ 5 — lOh ^ kpc. 

As an example, we again consider halo Aq-D-2 and show 
in the left panel of Fig. [2] the Y-projection of the surface 
mass fraction in subhaloes within r2oo. The right panel in 
the same figure shows the surface mass fraction of subhaloes 
averaged within azimuthal annuli as a function of the nor- 
malised radius R/r2oo- It is clear that the scatter in the 
projected mass fraction of subhaloes among different haloes 
is large. Within 0.1 r2oo, the mean fraction is ~ 0.005, with a 
scatter of a fact or of 10. The red l ine in the same panel shows 
the results from lMao et al.l (|2004l ). which were obtained from 
12 haloes (of galactic, group and cluster masses) and 30 
random projections. Their result lies somewhat higher than 
found here although still within the scatter. This is proba- 
bly due to the inclusion of group- and cluster-sized haloes 
in the averaging, which tend to have a higher substructure 
fraction due to their later formation times. The blue point 
indicates the required substructure mass fraction found by 
iDalal fc Kocha nck (2002) to be 0.02 (median, ranging from 
0.006 to 0.07 at 90% confidence). 



2.2 Adding "light" to dark matter haloes 

We put the source redshift Zs at 3.0. This is reasonable since 
many lensed quasars are at similar redshift. The lensing crit- 
ical surface density is given by 



(2) 



where Da, Dd, and Dda are the angular diameter distances 
between the source and the observer, the lens and the ob- 
server, and the source and the lens, respectively. For our 



adopted source and lens redshifts, Ecr = 1.82 x 10^ Mq 
kpc-^ = 7.95 X 10^° Mq arcsec-^ 

To produce multiple images, the maximum surface den- 
sity of a halo usually has to be super-critical. The left panel 
of Fig. [3] shows the surface density distribution for the halo 
Aq-D-2 projected along the Y-axis. Clearly, the central sur- 
face density of the (initial) NFW dark matter halo is below 
the critical value, and thus generally no multiple images ca n 
be produced (e.g. IWilliams et al.|[T999l : iRusin fc Mall200ll ) . 
This is hardly surprising, since for galaxy-scale strong lens- 
ing, the images form only a few kpc (projected) from the 
centre where baryons play a crucial role. Thus, one must 
incorporate a baryonic component in order to model the 
lensing galaxies more realistically, a topic we turn to next. 

Most gravitational lenses are early-type (elliptical) 
galaxies rather than late-type (disk) galaxies, as the for- 
mer are more massiv e and dominate the lensing cross- 
sections ([Turner et al.l [l984l ). There have b een many hy- 
brid models used for the l ensing galaxies (e.g. 'Keeton"200l|; 
Kochanck & White, I2OOII : [Oguri 2002; Jiang & Kochaney 
I2OO7I ). We use the spherical Hernquist profile to model the 
light distribution, since it approximates the de Vaucouleur's 
profile that has been observed for elliptical galaxies and 
bulges, and it has many known, convenient analytical prop- 
erties. 

The three-dimensional density and mass profiles p_ff(r), 
MH(r) for the Hernquist distribution are given by 
(|Hernquistlll990l ): 



PH{r) = 



2nr (r + a)^ ' 



MH{<r)^ AL 



(3) 



(r + a)2 



where Af* is the total baryonic mass, and a is a scale length 
related to the effective spherical radius (within which half 
of the mass is contained) by a = Tc/ (\/2 -|- 1). 

The profile is specified by two parameters a (or re) and 
A4i,, which are linked with the dark matter halo parameters 
r2oo and AA200 by 



M200 = Mdm + M* 



(4) 



r2oo ' ' M200 ' 
Notice that the mass of the main halo dark matter A^dm is 
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Figure 2. The left panel shows a contour map of the subhalo surface mass density fraction, which is the ratio of the surface mass in 
subhaloes to that in the total halo, for Aq-D-2 projected along the Y-axis. The right panel shows the mean distribution of subhalo surface 
mass fraction as a function of R/r 200, averaged over the three independent projections of each of the six Aquarius haloes at r edshif t 
2 = 0. The error bars indicate the 68% scatter among different projections and haloes. The red lines show the fit fr o m [Mao et al] l|2004h . 
The blue point indicates the median and 90% confidence level of the required fraction found bv lDalal &: KochanetJ 1I2OO2I ) (assuming the 
Einstein radius to be 0.02 r20o)- 



reduced by a factor of (1-/*) to conserve the total mass and 
the mass of substructures. 

The inclusion of the baryonic component affects the dis- 
tribution of the dark matter halo. Many studies have shown 
that the adjustment of the dark mat ter halo can be approx- 
imated by an adiabatic contraction (IBarnes fc Whitelll984l : 
iBlumenthal et al.lll986h . lGnedin et al. I l|2004h have proposed 
a modification to this simple model in order to take into ac- 
count the fact that particle orbits in realistic halos are not 
circular, but it is not clear whether this modification is able 
to repro duce accurately th e results of numerical simulations 
(see, e.g. lAbadi et al.ll2009l '). In view of this, we ha ve decided 
tofoUow, for simplicity, the procedure outlined bv lMo et al.l 
l|l998ll . Assuming that both the baryon and dark matter 
components follow an NFW distribution initially, baryons 
(/* percent of the total matter) then cool to form the galaxy 
at the centre, which causes the dark matter halo to contract 
adiabatically. After the adiabatic contraction, the dark halo 
follows a new profile and hosts a Hernquist galaxy at its 
centre. Note that we contract all the particles in different 
components (i.e., diffuse dark matter and subhaloes) in the 
same way. 

The two parameters (/ro and /*) are chosen accord- 
ing to two criteria (after adding the baryonic galaxy 
and accounting for the adiabatic contraction): (1) the 
projected dark-matter mass fractions inside the Einstein 
radii of the host galaxy haloes should range from 0.4 
- 0.7 l|Treu fc KoopmansI |2004 ): (2) the projected slopes 
are close to i sothermal at a few kpc from the galac- 
tic centre (e.g. Rusin et al] 20031: iRusin fc Kochanekll200l 



iKoopmans et al" I I2OO6I : iGavazzi et al.|[2"007h . or equivalently, 
the final rotation curves are roughly fiat from a few kpc out 



to a few tens of kpc (see Fig. |3J . Furthermore, /* should be 
smaller than the uni versal baryonic fract ion of ~ 17.5 per 
cent (from WMAP-3. [Spergel et al.ll2007h . 

We find that /ro = 0.05 and /* =0.1 satisfy these cri- 
teria well. From Fig. [3] (the left panel), it is clearly seen 
that after inserting the baryonic galaxy and taking the adi- 
abatic contraction into account, the total surface density is 
now super-critical and the corresponding Einstein radius is 
of the order of a few kpc, similar to that in many gravi- 
tational lenses. Notice however that our procedure is not 
self-consistent dynamically, since the inclusion of a bary- 
onic component will affect the evolution and survival of sub- 
haloes. We shall return to this point briefly in the discussion. 



3 LENSING METHODOLOGY 

Ai'-body simulations provide us with the positions (and ve- 
locities) of particles. For lensing calculations, we first project 
the particles onto a mesh in the lens plane (and tabulate the 
stellar surface density, and then smooth the surface density 
field appropriately. Using the smoothed surface density map, 
we can numerically calculate the lensing potential, defiection 
angles and magnifications. The details of the numerical pro- 
cedure are given in ii3.ll 

We test the accuracy of our numerical procedure by 
comparing with known analytical results, using a singular 
isothermal sphere realised through Monte Carlo simulations 
in ii3.2l We then relax the spherical assumption, and further 
test our procedure with an isothermal ellipsoid generated 
with a similar number of particles as those in the Aquarius 
simulations; the comparison results are presented in ii3.3l 
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Figure 3. The halo Aq-D-2: the left panel shows the surface density profiles S(-R) projected in the y-direction, and normalised to the 
critical surface density. Profiles are for cases before and after adding a Hernquist galaxy and the dark matter halo's adiabatic contraction, 
assuming the added baryonic component has 10% of the total mass and an effective radius of 5% of the halo virial radius (/* = 0.1, 
fic = 0.05). Line symbols are labelled inside the figure. The isothermal slope oc R~^) is indicated by the red line at the top right 

(see ^2.2l for details). The middle panel shows the mass distributions M{< r). The right panel shows the rotation curves Vc{r). The final 
total rotation curve is fiat from ~ 5h~^ kpc out to a few tens of kpc. 



3.1 From particles to lensing images 

3.1.1 Coarse and fine particle meshes 

We use a Particle-Mesh (PM) code for the lensing poten- 
tial calculation. The application of Fast Fourier Transforms 
(FFT) in the PM algorithm makes it computationally effi- 
cient. However it is limited in resolution by the finite mesh 
size and so cannot accurately represent regions with rapid 
density variations on the scale of the grid size. To increase 
the accuracy in the regions of interest (within a few kpc from 
the centres of galaxies), we establish two two-dimensional 
[2D] meshes: a coarse grid used for the potential field gener- 
ated by the mass projected outside the central {20h~^ kpc)^ 
region, and a fine grid for the mass within. Both grids have 
1024 X 1024 pixels, covering (4r2oo)^ and (40/i~^ kpc)^ (see 
§3.1.3) with resolutions ~ 0.6h^^ kpc and 0.04:h~^ kpc for 
the coarse and fine grids, respectively (the factor of two in- 
crease in the box size is due to the isolated boundary con- 
dition, see i33.1.3p . This resolution ensures that the tangen- 
tial critical curves are resolved with sufficient accuracy. In 
contrast, the inner radial critical curves may not be well re- 
produced, due to the finite resolution of the mesh. However, 
this is not a major concern since all the bright images that 
we are interested in form close to the outer (tangential) crit- 
ical curves. Furthermore, the resolution of the fine mesh is 
similar to the softening length of the simulations, and the 
density distributions in the very central regions are not ac- 
curately modelled in the simulations on smaller scales than 
the gravitational softening in the first place. 



3.1.2 Particle assignment with smoothed particle 
hydrodynamics kernel 

The surface density maps of the Aquarius haloes are ob- 
tained by assigning particles to the potential meshes us- 
ing the Smooth Particle Hydrodynamics (SPH) kernel 
l|Monaghanlll992l ). Although, in the end, we will approx- 
imate the underlying mass distributions of the Aquarius 
haloes by isothermal ellipsoids in order to circumvent prob- 



lems caused by discreteness noise (see H3.3p . SPH-smoothed 
density fields are used as an intermediate step to generate 
basic lensing properties (e.g. critical curves and caustics) to 
constrain the best-fit isothermal ellipsoids. For more detail, 
see !j4l 

The advantage of the SPH assignment is that it adjusts 
the smoothing scale according to the local density environ- 
ment: particles in a high density region are mildly smoothed 
while those in a low density region are smoothed more. For 
each particle, a smoothing length Ti is calculated according 
to the local number density in its 3D neighbourhood. The 
particle mass is then assigned to all the mesh cells that are 
within a circle of 2TC in radius in its neighbourhood. The 
3D density kernel can be integrated along the line of sight 
analytically to obtain the surface density distribution: 



^[-(8 + 52u'^)Vl-u'2 + (16 + 26m^)^ 
-Qm" lnu + 3^2(16 + 4^2) ln(l 
-311^(16 + 
i[2V4^(l 



vr 



ln(2 + V4-m2)], ifl>tt>0 
01, if2>it>l 



-3it2(16 + u2)ln(2 + 
0, if M > 2 

(5) 

where u = r/Ti is the distance from the cell centre to the 
particle normalised to Tl, and the total mass within u < 2 
is unity. 

The smoothing length Ti for each particle depends on 
its local density and is controlled by the parameter A^'ngb, 
the number of particles that are contained within radius 
TC. A good smoothing procedure should reduce the numeri- 
cal noise without smoothing excessively out the real density 
fiuctuations (e.g. substructures). 

The total mass that has been assigned to the neigh- 
bouring cells should be equal to the mass of the particle. 
However this is only approximately true due to the discrete- 
ness of cells. In particular, mass conservation is quite poorly 
observed when the smoothing length TC is only a few mesh 
cells, which may happen in a dense environment. For par- 
ticles with 2Tt < 15 cell sizes (30 cells in diameter), we 
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therefore renormalise each individual kernel so that the to- 
tal mass is conserved during the assignment. 

We find in practice that SPH assignment is superior to 
Cloud-In-Cell (CIC) assignment in terms of reducing dis- 
creteness noise. For a singular isothermal sphere realised 
with 10*^ particles, the SPH-smoothed (iVngb = 32) and CIC- 
smoothed surface density fields show fluctuations of 2% and 
30% relative to the analytical results, respectively. For a re- 
alisation with 10^ particles, the fiuctuations decrease to 1% 
for the SPH assignment (A'ngb = 320, with the same smooth- 
ing length, n oc N~'^'^nI'^1 \U et al.ll200'^ ) and to 10% for 
the CIC assignment. 



3.1.3 Isolated boundary conditions 

Periodic boundary conditions are most natural for Fourier 
Trans forms, but are not appropriate for lensing galaxies. We 
follow iHocknev &: Eastwood! l|l98ll ) to eliminate the (alias- 
ing) effects due to "mirror" particles by using a mesh twice 
as big as the simulated lens system, padding the region out- 
side the simulation volume with zeros. A truncated 2D grav- 
itational force kernel is tabulated onto the same simulation 
mesh, and then convolved with the assigned surface den- 
sity field. The gravitational effect is accurately reproduced 
within the region where the mass has been distributed (See 
iHocknev fc Eastwood! [l98l] . for more technical details). We 
adopt this procedure throughout this work. 



3.1.^ Lensing potential, deflection angle and magnification 

After the discretisation of the surface density field through 
SPH assignment and the tabulation of the truncated 2D 
gravitational kernels on the meshes, the potentials and their 
derivatives are easily calculated by convolutions which can 
be efficiently implemented in Fourier space. 

In particular, the effective lensing potential ^i){6) is the 
convolution of the surface density T,{9) and the 2D kernel 
Inlel: 



T.{e')\u\9-e'\d^e'. 



(6) 



The deflection angle Q.{6) is the first derivative of the lensing 
potential, V'(^)) ^^nd is thus the convolution of the surface 
density E(e) and the 2D force kernel e/\e\^: 



E(e') 



■ d' 



(7) 



The convergence k{9) (the surface density normalised to Ecr) 
and the shear "/{6) are second-order derivatives of the lensing 
potential tpiO): 

«: = (V'll + ^22)/2, 71 = (l/"!! - V'22)/2, 



72 = V'12 = V'21, 7^ = 7l + 72, 



(8) 



where the derivatives are taken with respective to the index 
1 (x) and 2 (y). Numerically, the convergence and shear can 
be calculated through 4th-order finite differencing from the 
defiection angle a{9). The magnification /i(S) is related to 
the convergence and shear by 



3.1.5 Image finding and cusp relation 

Since all the lensing quantities are now known, it is straight- 
forward to find the images for any given source position. To 
this end, we construct a separate mesh in the image plane, 
with a resolution (0.02/i~^ kpc) higher than the fine po- 
tential mesh discussed in H3.1.1I the lensing properties (de- 
flection angle, magnification etc.) on this ultra-fine mesh are 
found through bi-linear interpolation. We then search image 
positions (and magnificatio ns) using the N ewton-Raphson 
and triangulation methods l|Schneider et al.i.l992 ). 

Of particular interest to gravitational lensing are the 
critical curves and caustics. Critical curves in the image 
plane are a set of points where the magnification is formally 
infinite for a point source, ^(9) — !• oo. In practice, they are 
identified according to the fact that the magnifications have 
different signs (i.e., different parities) for images on differ- 
ent sides of a critical curve. Critical curves are mapped into 
caustics in the source plane, which can be easily obtained 
through the lens equation. 

Most strong lenses occur in elliptical galaxies since 
they have larger len sing cross-sections than spiral galaxies 
(|Turner et al.![l984l ). They typically form two distinct sets 
of critical curves and corresponding caustics: the tangen- 
tial ("outer") and radial ("inner") critical curves, which are 
mapped into tangential ( "inner" ) and radial ( "outer" ) caus- 
tics (see Fig. (51 for an example). A source inside the central 
caustic usually produces five images: four close to the tan- 
gential critical curve and one central image which is usually 
too faint to be observable (and is of no interest to us for the 
present work). 

We are particularly interested in sources that are close 
to the cusps of the tangential caustic ("cusp sources"). For 
cusp sources, three close images form around the tangential 
critical line, with alternate parities. There are two different 
kinds of cusp sources and corresponding image configura- 
tions. As illustrated in Fig. (5] a "major cusp" source forms 
three images around the tangential critical curve on the same 
side of the source (with respect to the centre of the lens) 
while a "minor cusp" source forms three close images on the 
opposite side of the source. 

In any smooth lensing potential, for a source very 
close to a cusp, the three close images satisfy an 
asym ptotic magnification relation (the "cusp-caustic rela- 
tion" ; !Blandford fc Naravan![l986l: !Schneider fc Weiss!!l99A 



!Keeton et a"i]!2003! ): 



-Rcusp 



0, 



(10) 



with the total absolute magnification |/^a| + |ms| + Ipc| —* oo- 
For each of the cusp sources, we define an image open- 
ing angle A^, ranging from to tt, which measures the angle 
(from the lens centre) of the outer images of the close triple. 
Notice that both A6 and i?cusp are observable. In a smooth 
lens potential, as a source moves to a cusp caustic, both /S.9 
and -Rcusp decrease asymptotically to zero. As can be seen 
from Fig. (6] there are two leading patterns on the i?cusp — AS 
diagram due to "major" and "minor" cusp sources. Gener- 
ally speaking, the major cusp sources have larger iicusp than 
the minor cusp sources for the same image opening angle. 

The cusp-caustic relation predicts that in smooth lens 
models -Rcusp would asymptotically approach zero when a 
source moves towards the caustic. However, the presence 
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of (clumpy) substructures will break down the smooth po- 
tential assumption in the asymptotic cusp-caustic relation, 
resulting in substantial deviations in -Rcusp values and other 
quantities (such as image positions and time delays) from 
simple predictions. Therefore the examination of the cusp- 
caustic relation is a way to test for the presence of sub- 
structures that are projected near the (tangential) critical 
curves. However, caution must be exercised because, even 
for smooth lens models, a high iicusp is possible. There are 
many factors that affect the -Rcusp distribution apart from 
the presence of substructures, e.g. the mass distribution of 
the lens (radial profile and the ellipticity) , external shear 
from the environment, and the s election criteria of t he cusp 
sources (for more discussions see iKeeton et allbooal ). 

3.2 Singular isothermal sphere 

We test our lensing simulation code with Monte-Carlo real- 
isations of singular isothermal spheres (SIS), for which an- 
alytical lensing properties are known. Each of our SIS con- 
tains a mass of 10^^ Mq within a virial radius of WOh~^ 
kpc, realised with 10® and 10* particles; the SPH assignment 
parameter is chosen to be A'^ngb = 32 and A'ngb = 640 for 
the two cases respectively. Fig. U shows the numerical ac- 
curacy of the deflection angle, convergence (surface density) 
and magnification in our numerical procedures. For the 10® 
particle case, two Monte-Carlo realisations (cyan and blue 
curves) are shown. For the 10* particle realisation, the un- 
certainties around the Einstein radius (at about 0.02 r2oo, 
defined by a total magnification fi{9) > 20) are 0.2% for 
the defiection angle, 1% for the convergence, and < 10% for 
the lensing magnification (estimated by the half width half 
maximum of the probability distributions). The deviation 
towards the centre is due to the fact that the finite mesh 
resolution of the Particle-Mesh code fails to represent the 
singular behaviour at the centre of the SIS. The significant 
deviation of the magnification seen near the Einstein radius 
is due to the divergent behaviour of the magnification close 
to the critical curve, where fj, — 1/((1 — k)^ ~7^) — * 
when AC = 7 = 0.5 at the Einstein radius for the singular 
isothermal sphere. 

3.3 High-resolution isothermal ellipsoid 

We simulate an isothermal ellipsoid (IE) with 10® and 10* 
particles (as in the Aquarius haloes) . Such an isothermal el- 
lipsoidal distribution is modelled as an oblate spheroid with 
axis ratio gs and with a density distribution: 

p(x(s'o+R' + zV<ilr\ (11) 

where So is a core radius, and {R, z) are the cylindri- 
cal coordinates. It i s spe cified by three parameters (see 
iKeeton fc Kochanekl 1 19981 for details): the effective criti- 
cal radius bi, the eccentricity of the mass distribution e = 
(1 ~ QzY^^ ^ Build a core radius So. The parameters for the 
isothermal ellipsoidal halo are adjusted so that its criti- 
cal curves and caustics match those for the halo Aq-F-2 in 
the z-projection. The parameters are bi = 0.4", — 0.8, 
So = 0.1" and the major-axis of the surface density ellipse 
is rotated by ~ 7r/8 with respect to the X-axis. 

Fig. [5] shows analytical and numerical critical curves. 



caustics and -Rcusp maps for sources located inside the dia- 
mond caustics. The two critical curves nearly overlap each 
other, with barely noticeable wiggles in the numerical result. 
The caustics also agree reasonably well for the 10*-particle 
case, but for the 10®-particle realisation (blue lines), higher- 
order singularities such as swallowtails are clearly seen close 
to the cusps and along the fold caustics. These arise due 
to numerical noise. The numerical -Rcusp map shows visible 
distortions compared with the smooth contours in the ana- 
lytical results, even in the central region. 

Fig. |6] presents the -Rcusp distributions for caustic 
sources (indicated in the left panel) with image opening an- 
gle A9 < 90° . The analytical results show two distinct peaks 
due to major-cusp and minor-cusp sources with a sharp 
dropoff around -Rcusp ^ 0.12. In contrast, for the numeri- 
cal distributions, even the 10*-particle realisation shows a 
much broader profile than the analytical one, with an ex- 
tended tail out to -Rcusp ~ 0.4 due to numerical noise. No- 
tice that the numerical noise behaves in a similar way as real 
substructures in the -Rcusp distribution. 

In the current numerical set-up with -/Vngb ~ 640 for a 
10*-particle simulation, the smoothing length Ti. in the cen- 
tral region (of the halo) approximately reaches the softening 
length of the Aquarius simulation, also roughly the cell res- 
olution of the fine mesh. Increasing TL will indeed further 
suppress the noise, but it may also over-smooth the under- 
lying density field. Below we outline an alternative way to 
approximate the smooth underlying density fields for the 
host galaxy haloes. 



4 RESULTS 

4.1 Lensing Predictions for Aquarius Haloes 

In this section, we will apply our lensing methodology to 
the Aquarius haloes (and a baryonic component modelled 
as a Hernquist profile), and study the violation of the 
cusp-caustic relation due to the dark matter substructures 
therein. In principle, we should compare the lensing proper- 
ties of the simulated haloes with and without substructures 
in order to assess the effects of substructure. However, Fig. 
[6] shows a substantial broadening of the -Rcusp distribution 
due to numerical noise, which will significantly confuse sig- 
nals from real substructures. 

As mentioned above, the total (dark matter plus 
baryons) mass profile of each halo is adjusted to resemble an 
isothermal distribution in the central region. To avoid exces- 
sive discreteness noise, we go one step further and adopt a 
fitted isothermal ellipsoid (with a density distribution given 
by eq. [TT]) rather than the original particle distribution 
for subsequent lensing calculations. The parameters of the 
isothermal ellipsoid model are adjusted to match the orig- 
inal critical curves and the caustics of each Aquarius halo 
(together with the central galaxy) in a particular projec- 
tion. We add the particle distributions of substructure in the 
Aquarius simulations (assigned to meshes according to the 
CIC algorithm) to the isothermal ellipsoid that fitted to the 
main galaxy halo, and compare its lensing properties with 
those of the smooth underlying isothermal ellipsoid. In this 
approach, the analytical solutions of the fitted isothermal 
ellipsoidal potential and its derivatives are tabulated on the 
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Figure 4. The numerical accuracy of the deflection angle, the convergence, and the magnification for Monte-Carlo realisations of singular 
isothermal spheres. The top panels show the ratios of the numerical to the analytical results as a function of radius. The deviation in 
the numerical magnification (on the right) towards the centre is due to the finite mesh resolution of the Particle-Mesh code, and that 
seen near the Einstein radius (at about 0.02 r20o) is due to the divergent behaviour of the magnification close to the critical curve. 
The corresponding probability distributions for images with a total magnification above 20 (around r ~ 0.02 r2oo) S'l's presented in the 
bottom panels. The cyan and blue curves are for two 10^-particle realisations (with W^gb = 32) while the red curve is for a 10*-particle 
realisation (with A''ngb = 640). 



image grid, including the cusp-caustic relation. By doing so, 
no Poisson noise of the underlying main halo is introduced, 
and so any confusion to the results from substructures is 
avoided. 

We fit an isothermal ellipsoidal model to each of the 
three independent projections of each galaxy halo. The six 
fitting parameters are: (1) the efi'ective critical radius bi, 
(2) the axis-ratio of the surface density ellipse qa, (3) the 
core radius So, (4)-(5) the X- and y-offsets of the projected 
centre X^, Yc, and (6) the rotation angle of the major axis of 
the surface density ellipse RA. The uncertainties of the fitted 
parameters are Afoi = 0.003", Ags = 0.01, A5o = 0.001, 
AXc = AFc = 0.002", /\RA = 0.0047r. The relative errors in 
the fitted critical curves and caustics are < 10 % for diflerent 
projections of all the simulated galaxy haloes. 

In Table O we list the isothermal ellipsoid parameters 
of the main haloes (fei, 53 and So). The critical radius bi is 
of the order of 0.3" to 0.9". The separation between images 
(~ 2bi) is in the range of the obs erved gravitational lenses 
(which peaks around 1", see e.g. iBrowne et al.ll2003D . The 
axial ratios also match the observed lenses quite well. There 
is one exception, however. The core radius So is quite large, 
of the order of (0.05" — 0.1", corresponding to a few hundreds 
of pc). Such a core size is larger than those inferred from 
gravitational len ses which are in general consis t ent with zero 
core radius (e.g. IWallington fc Naravanlll993l : iRusin fc Mai 
l200ll : IOguri et al.ll200ll : lLi fc Qstrikej|2003l '). This is a direct 



result of the implementation of the Hernquist profiles, which 
follow a logarithmic density slope of — 1 in the central regions 
(see Fig. [ij . However, this artifact should have little effects 
on the images we are interested in, which are close to the 
outer critical curve. 

To examine the violation of the cusp-caustic relation, 
we generate about 10000 cusp sources in each case, and 
calculate the resulting 7?cusp distributions. All these cusp 
sources are inside the central caustic and close to the cusps, 
where the corresponding triple images have opening angles 
A9 < 90°. The results of all 7 studied Aquarius haloes (in 
21 projections) are given in Fig. [7] to Fig. 1131 As mentioned 
above, the probability density distribution of -Rcusp often 
shows two peaks for the smooth haloes which are produced 
by the major and minor cusps, respectively. For the "naked" 
cusp cases (where the central diamond caustic protrudes the 
outer elliptical caustic), the distributions of iicusp vs. the 
opening angle A9 are somewhat truncated below certain 
opening angles (see Fig.|5]for the halo Aq-C-2's y-projection 
for an example). Empirically, it is rare for massive lensing 
galaxies to prod uce naked cusps. There is only one candi- 
date APM08279 l|Lewis et al.ll2002l ). and that is hkely due to 
lensing by an edge-on spiral rather than an elliptical. The 
four naked cusp cases from our simulations are caused by 
the large cores in the central density profiles of the lensing 
galaxies. We exclude these four naked cusp cases in the final 
statistic calculations (their inclusion does not significantly 
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Figure 5. Critical curves, caustics, and cusp-caustic relation -Rcusp maps for a Monte-Carlo realisations of isothermal ellipsoid with 
A^p = 1.6 X 10** and Af^gb = 640. The top panels show the critical curves and the caustics. The position and corresponding images are 
shown for a "major cusp" (solid squares) and a "minor cusp" source (open squares). The bottom panels show the i?cusp maps from the 
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with image opening angles A9 < 90°. The left panel shows the source positions with respect to the tangential caustic. The middle panel 
shows both the numerical (10®-particle realisation) and analytical results of -Rcusp vs. A9. The right panel shows the probability density 
distributions of iJcusp for the analytical IE (red), and the two Monte-Carlo realised haloes with 10® (cyan) and 10* particles (blue), 
respectively. 
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alter our results). Strong violations of the cusp-caustic re- 
lation due to substructures are seen in some cases, e.g. for 
the y-projection of the Aq-B-2 halo (see Fig. |8]). However, 
most of these cases have small cusp-lensing cross-sections 
(listed in Table [3l Column 8), defined as the areas covered 
by cusp sources whose images satisfy AS < 90°. The mean 
probability of cusp violations calculated below are weighted 
by the cross-sections (see ea. ll2|) . As can be seen from Fig. [7] 
to Fig. 1121 the scatter in the cusp violation is large between 
different projections of different haloes. Also notice that the 
halo Aq-A-2 aX z = 0.6 (Fig. I13|l does not show a significant 
difference from the redshift-zero haloes in the violation of 
the cusp-caustic relation. 

To see which massive substructures cause the cusp- 
caustic violation, we calculate the -Rcusp distribution due 
to subhaloes more massive than Kf'h-^MQ, W'^H-^Mq, 
lO''/i"^M0, and lO^/i'^M©, respectively. Fig. [14] shows 
one typical example, for the halo Aq-D-2 along the Z- 
projection. We find that in most cases substructures with 
masses m-sub < 10^ to lVph~^ Mq dominate the contribution 
to the violations of the cusp-caustic relation (see the Col (9): 
Afsub.cr in Table [Sjl. Notice that previous studies on cusp vi- 
olations typically resolve haloes larger than ~ 10*/i~^Mq, 
and thus would not have been able to evaluate the effects of 
substructure accurately. However, the addition of subhaloes 
with rrisub ^ IQ^ Mq does not appear to increase the vio- 
lation frequency significantly (compare the three right pan- 
els). We return to the convergence issue as a function of 
subhalo mass in ^ 

Notice that most subhaloes that are projected close 
to the critical curves are due to chance alignment. Fig. 1151 
shows the spherical halocentric distance distribution for the 
subhaloes that are within a projected distance of 0.05 r2oo 
(~ 2.5 Einstein radii). The fractions of subhaloes that are 
physically located within a spherical radius of 0.05 r2oo are 
15%, 18%, 15% and 0% for subhaloes more massive than 
IQ^K'^Mq, 10^/i"^Mq, IQ'h-^MQ and 10*/i"^AfQ, respec- 
tively. The large median halocentric distances, ~ 0.2 r2oo in 
all cases, also show that projection effects are substantial. 



4.2 Comparison with observations 

iKeeton et all (|2003l ) summarised 19 published quadruply 
imaged systems. Seven of them are detected at radio wave- 
lengthfl' Radio lenses are free from dust extinction. Due to 
their large emission regions, they are less likely to be af- 
fected by microlensing. In contrast, microlensing is likely to 
affect optical/IR flux ratios and so we treat them differently 

below^ 

iDalal fc Kochanekl (|2002D studied seven four-image 
radio -lensing systen is: MG0414 -I- 0534 l|Hewitt et al.1 
Il992l). B0712-H472 jjackson et al.l Il998l). PG1115-F080 
JWeymann et al.lll98d). B1422-I-231 (iPatnaik fc Narasimhal 
l200lh. B1608-H65 6 l|Fassnacht et all Il996l). B1933-H503 



Table 2. The image opening angle and -Rcusp for the observed 
cusp-caustic lenses, taken from lAmara et aL I l l2006fl . 



Lens 


AS 


-f^cusp 


Band 


B0712-f-472 


79.8° 


0.26 ± 0.02 


radio 


B2045-I-265 


35.3° 


0.501 ± 0.035 


radio 


B1422-I-231 


74.9° 


0.187 ± 0.006 


radio 


RXJ1131 - 1231 


69.0° 


0.355 ± 0.015 


optical/IR 


RXJ0911 + 0551 


69.6° 


0.192 ± 0.011 


optical/IR 



and found that six show anomalous flux ratios, which might 
be due to the effects of substructure lensing. Among all 
the detected radio lenses, three (B0712-I-472, B1422-I-231 
and B2045-I-265) show a typical cusp-caustic geometry 
(with AS < 90°) and violations of the cusp-caustic relation. 
Another two lensing systems observed in the optical/IR 
band are also cusp-caust ic lenses with AS < 90° : RXJ1131- 
1231 llSluse et al.1 l2003h and RXJ0911-h0551 jBade et al.1 
1 19971 : iBurud et al.l 1 19981 ). Both have unexpected large 



by microlensing l|Morgan et al.ll2006l : lAnguita et aLlbOOSl ). 
Table [2] lists the -Rcusp and AS values for the five observed 
cusp-caustic lenses. Three out of the five cusp lenses are 
detected at radio wavelengths, thus their large -Rcusp values 
are unlikely due to microlensing. We treat these three radio 
lenses as cusp-caustic violations due to substructure lensing. 
Below we will calculate the probability for the simulations 
to reproduce such an observed violation rate. 

For each galaxy and each projection, we calculate the 
violation probability that the predicted 7?cusp is larger than 
the observed -Rcusp value 0.187 for B1422-I-231, which shows 
the smallest violation (smallest -Rcusp value) among the five 
cusp lenses with AS < 90°. The cross-section weighted vio- 
lation probability is given by 

Pa =^/a,.P,(i?cusp > 0.1871AS < 90°),/,,,, 



(12) 

where the summation i = 1, (21— 4) is for the seven haloes 
along the three independent projections of each, excluding 
the four naked cusp cases, and ai is the cross-section in the 
source plane for producing three close images with opening 
angle AS < 90° . Using the above formula, we flnd the mean 
probability pa ~ 6.4% for -Rcusp > 0.187. Notice that this 
probability estimate is only approxima te, since we have n ot 
considered the magnification bias fe.g. [Turner et al]|l984l ). 

To have three (radio) lensing cases with 7?cusp > 0.187 
(due to substructure lensing rather than microlensing) out 
of the five cusp lenses (AS < 90°) observed so far, the prob- 



ability is Cipt{l-pa 



2.3 X 10"''. The low probability 



( Svkes et al] 119981 ) and B2045+265 jFassnacht et al.lll999l ) 



suggests that the subhalo populations in the inner regions of 
the Aquarius haloes with Hernquist galaxies are insufficient 
to explain the observed frequency of flux anomalies in the 
cusp lenses. 



* B0 1 28-I-437 llPhillips et aLllioOol) . B0712-t -472 jjackson et al.l 
19981: iJackson et al.l | 2000|). B1422+231 lllmpey et al.l Il996t 
Patnaik fc Narasimha' '200l'). B155 5-I-375 l lMarlow et al l 
19991). B16 08+656 (Koopmans fc Fass nachd Il999l). B193 3-I-503 
llCohn et al.ll200ll) and B2045+265 iPassnacIit et aI.lll999^ . 



5 DISCUSSION AND CONCLUSIONS 

In this paper, we have used the ultra-high resolution Aquar- 
ius simulations to study the effects of substructure lens- 
ing. We incorporate the effects of baryons in the main 
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halo by adding a stellar component (modelled as a Hern- 
quist profile), and then take into account its eflFects on the 
dark matter halo through adiabatic contraction. The den- 
sity profiles and lensing properties except the flux ratios are 
broadly consistent with the observed gravitational lenses. 
Using Monte Carlo simulations, we find large numerical 
noise for an isothermal halo populated with 10* particles, 
which shows considerable scatter in the i?cusp distribution 
for cusp lenses. In the end, we therefore study the substruc- 
ture lensing by modelling the smooth underlying galaxy halo 
as an isothermal ellipsoid and superimposing the subhalo 
population from the Aquarius simulations. In this way, we 
focus on the lensing effects of subhaloes and avoid any con- 
fusion from numerical noise in the A*'-body realisation of the 
simulated main haloes. 

Our study finds that even with the much better re- 
solved subhalo population of the Aquarius simulations, the 
observed cusp lenses still violate the cusp-caustic relation 
more frequently than predicted by A'^-body simulations. 

The Aquarius haloes are Milky Way type haloes in 
terms of their masses, while many lenses are ellipticals, 
which are more massive. Among the five cusp lenses we com- 
pare our resuhs with, three of them (B2045-I-265, RXJ1131- 
1231, RXJ0911) are more massive than our simulated haloes 
and have Einstein radii twice as large as those of our haloes. 
The other two lenses (B0712-I-472, B1422-I-231) have Ein- 
stein radii and velocities (circular velocity or velocity dis- 
persion) roughly comparable to the relatively massive haloes 
in the Aquarius simulations. As shown in Fig. [2] (the right 
panel) the projected subhalo mass fraction increases with 
the projected radius R. If we have under-estimated the Ein- 
stein radii bj (e.g. because of uncertainties in the addition 
of the central galaxies), we could have potentially under- 
estimated the violation rates due to the lack of enough sub- 
structures at smaller radii. We artificially increase the Ein- 
stein radii of the simulated haloes by a factor of two to study 
the violation probabilities due to a higher fraction of sub- 
structures at larger radii. The mean subhalo mass fraction 
within a 0.l"-annulus around the new Einstein radius would 
increase from /sub.annu ~ 0.19% to 0.24%, and the mean vio- 
lation probability would increase from ~ 6.4% to 14.0%. 
The probability of reproducing the observed violation rate 
would increase from 0.2% to 2%. 

Another concern is that due to the finite particle mass in 
A''-body simulations, the central cusps of the subhaloes may 
not be resolved, which may potentially result in an under- 
estimation of the ability of the subhaloes to induce pertur- 
bations to the lensing potential. We consider an extreme 
case assuming all subhaloes are point-like sources with their 
masses and locations from the simulations. In this scenario, 
/sub.annu roughly remains at 0.18%, however, po- increases to 
15.1%. The probability to reproduce the observed violation 
rate increases to 2.5%. 

These low probabilities suggest that the subhalo popu- 
lations in the central regions of the Aquarius haloes are not 
sufficient to explain the observed frequency of violations of 
the cusp-caustic relations. It is important to ask whether our 
results will change significantly if even lower-mass subhaloes 
are resolved. We argue that this is unlikely to be the case. 
The total subhalo lensing cross-section is an integral of the 
cross-section of subhaloes of each mass weighted by their 
abundance. As shown in ^ most of the perturbing sub- 



haloes have relatively low mass (msub < 10^ to 10^h-'^M@). 
Their abundance scales as dA^(msub)/dmsub oc rn'^^^^ . For a 
galaxy (subhalo) approximated b y a SIS, the l ensing cross- 
section roughly scales as ct* (e.g. Turner et al.|[l984l ) where 
a is the one-dimensional velocity dispersion. For Aquarius 
subhaloes, mgub oc V^ax l|Springel et al.l [ioosi ') . where Vmax 
is the maximum circular velocity. If cr oc Vmax, then the in- 
tegrated lensing cross-section will be oc m^[^. On the other 
hand, for a point lens or an elliptical galaxy, the lensing 
cross-section is proportional to the lens mass, and the inte- 
grated lensing cross-section would be oc m'^;^^. In all these 
cases, the subhalo lensing cross-sections are biased towards 
relatively massive subhaloes in the projected central region, 
and the incorporation of even lower mass subhaloes should 
not change our results significantly. 

We mention in passing that a warm dark matter sce- 
nario would suppress the formation of small subhaloes, mak- 
ing it even m ore difficult to explain the observed cusp viola- 
tions (see e.g. lMiranda fc Macci6ll2007^ . Below, we compare 
our study with previous work, before discussing its limita- 
tions and outlining possible future work. 



5.1 Comparison with previous studies 

There have been a number of studies of substructure lensing 
using numerically simulated haloes, including those from hy- 
drodynamical simulations. Below we compare a few of these 
stu dies with our own. 

iDalal fc Kochanekl l|2002h concluded that at the 90% 
confidence level, a substructure fraction of 0.6% to 7% can 
explain the observed anomalous fiux ratio. For the Aquarius 
subhalo population. Table |3] Col (5): /aub,annu shows such 
fraction averaged over a thin annulus around the outer tan- 
gential curve, which is always below 1%, sometimes much 
smaller (not to be confused with /sub in Table [1] and Fig. [2j 
which refers to the subhalo mass fraction within r2oo)- This 
is the primary reason why our predicted cusp violations are 
sma ller than the observe d violation frequencies. 

iBradac et al.l ||2004 ) used hydrodynamical simulations 
of ISteinmetz fc Navarrol (| 2002^) and concluded that the pre- 
dicted cusp violations due to substructure are comparable 
to those observed. Their simulated halo has ~ 10^ parti- 
cles, resolving subhaloes down to 5 x 10* Mq. As the au- 
thors pointed out, the numerical noise may be as high as 
5%. The observed high-order singularities in their simula- 
tions are much higher than ours (comparable to the caustic 
structure shown in Fig.[5]for lO'^ particles). It is possible that 
their high numerical noise may have produced too many ar- 
tificial violations, although we note that they used Voroni 
density estimation to reduce the discreteness noise. 

Our conclusion that the dark matter subhalo popula- 
tion may be insufficient t o explain t h e obs e rved cusp vi- 
olations is consistent with Mao ot al. (2004"), 'Amara et ajj 
( 20061) . iMaccio et al.1 |2006) and Maccio fc M iranda (200d). 
The number of particles used in those studies is roughly 
two orders of magnitude smaller than here. In particular, 
the study by iMao et al] l|2004h found large scatter among 
different haloes, a conclusion confirmed by our results. 
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5.2 Limitations of the present study and future 
work 

The most severe limitation of our study is that the high- 
resolution simulations used here include only dark matter. 
Without baryons, these haloes are sub-critical (see ^ and 
incapable of producing multiple images. We are therefore 
forced to incorporate a model for the baryonic galaxy at the 
centre of each halo. The galaxy changes not only the over- 
all dark matter profiles (taken into account by adiabatic 
contraction), but also the dynamical evolution of subhaloes, 
an effect which is not considered here. On the one hand, 
the increased baryonic density at the centre of the halo will 
make the subhaloes feel stronger tidal forces, particularly 
those that come close to the centre. On the other hand, the 
baryons within subhaloes will make them more resilient to 
tidal disruption. It is not clear which effect will dominate. 
We comment, however, that the subhaloes that come very 
close to the centre may have already been tidally stripped or 
disrupted, and thus most of the surviving subhaloes that can 
be identified by SUBFIND may have quite large peri-centre 
passages. As a result, the effects of baryons in the host halo 
may not change the results very significantly. However, we 
caution that, SUBFIND, like most substructure finders, has 
difficulties in identifying subhaloes in the densest regions 
of the halo and assigning them correct masses. Empirically 
the Milky Way does not seem to host many luminous satel- 
lites close to the centre. Hydrodynamical simulations can 
in principle address this issue directly (subje ct to the un- 
certai nties in the treatment of gas processes). iMaccio et al.l 
found a factor of two increase in the number of sur- 
viving satellite galaxies (with masses above 10^ M©) in the 
centres of galaxies when including baryons in the CDM sim- 
ulations, but concluded that even this was not sufficient to 
explain the flux anomaly problem. 

Observationally, it is interesting that more than one 
half of the CLASS lenses a ppear to show luminous compan- 
ion g alaxies in projection l|Brvan et al.ll2008l : I Jackson et al.l 
12009 ^'). and their inclusion in the models appears to alleviate 
the anomalous flux ratio problem (see below). This may be 
just a statistical fluke due to the small sample size (22 lenses 
in total) or some of th ese may be due to chance alignment 
along the line-of-sii ^ht lichen et al.ll2003l: l\yambsganss et al.l 
l2005l : lMetcai3l2005l a.b: ]Miranda fc Macci6ll200'it l. Neverthe- 
less, for the three radio lenses that show apparent cusp vi- 
olations (see Table [5} , the most serious cas e is B2 045-f265 
with Ticuap ~ 0.5. Recently, iMcKean et all (|2007h found a 
galaxy, G2, which is about 0.66 " away from the main lens- 
ing galaxy Gl (at redshift 0.867), and about 3.6 to 4.5 mag- 
nitudes fainter than Gl depending on the wavelength. The 
photometric redshift of G2 is consistent with that of Gl (al- 
though also consistent with a redshift ~ 4—5). The inclusion 
of this faint satellite galaxy in the model can explain the flux 
anomaly reasonably well, although the satellite is required 
to be very flattened with an axis ratio of 8:1, which may 
not be realistic. This case highlights the potential roles that 
the luminous satellites may play in the anomalous flux ra- 
tio problfiim_We note, however, that numerical simulations 
bv lDolag et al.l l|2008l ) showed that star-dominated galaxies 
(not traced by dark matter only simulations) appear to con- 
tribute only ~ 10% of the subhalo population in clusters 
of galaxies. It is unclear however whether this cluster-based 



result can be extrapolated to galaxy scales where cooling is 
more efficient. We plan to use semi-analytical galaxy cata- 
logues in the Aquarius simulations to address this issue more 
quantitatively in subsequent work. 

Substructures not only perturb the flux ratios, but also 
affect the image positions. In Table [31 we show the maxi- 
mum perturbation of the deflection angle, aaub.max, within 
the central 2" x 2" region, produced by all the subhaloes 
within r2oo . The maximum deviations range from a few milli- 
arcseconds to < 0.1 arcseconds. They may leave observable 
signatures on clos e pair imag es such as that observe d in 
MG2016-f 112 (Ko opmans et a l. 2002; More et al.ll2009l ). We 
find that most of these astrometric deviations are dominated 
by one large, nearby subhalo. This clearly warrants further 
work in the near future. 
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Figure 7. Lensing properties for the halo Aq-A-2, in three independent projections. The top panels show the critical curves (red) and 
caustics (blue) superimposed on top of the subhalo population. The middle panels show Rcusp vs. the image opening angle A9. Large 
R cusp values (red) are due to substructure. The triangle pattern (green) gives predictions for the smooth counterparts. The bottom 
panels show the corresponding probability distribution functions (PDFs) of -Rcusp for cusp sources with AO < 90°. The violation of the 
cusp-caustic relation can be seen from the excess of iJcusp at large values (red) over the smooth counterpart curve (green). The iJcusp 
values for the three radio and two optical/IR cusp lenses are indicated by vertical solid and dashed bars (see Table [2]l. 
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Figure 8. For the halo Aq-B-2, the symbols arc the same as in Fig. [7] The truncated triangle pattern in the Z-projcction is due to 
naked cusps of the central caustic. The strong violation of the cusp-caustic relation seen in the Y-projection is caused by subhaloes with 
"isub < lO*/i-lM0 with a violation rate P(Kcusp > 0.187) = 64%. 
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Figure 9. For the halo Aq-C-2, the symbols are the same as in Fig.[T] The truncated triangle pattern in the y-projection is due to naked 
cusps of the central caustic. The halo in this projection has large ellipticity, which results in large -Rcusp values. The strong violation in 
the X-projection is caused by subhaloes with m^ub < I0^h~^ Mq with a violation rate P(-Rcuap > 0.187) = 19%. 
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Figure 10. For the halo Aq-D-2, the symbols are the same as in Fig. [7] The strong violations in the Y- and Z-projection are caused by 
subhaloes with rrisub < lO"^ Mq and m^^^, < 10^h~^ Mq, respectively, with violation rates P{Rcubp > 0.187) = 9.7% (y-projection) 
and 31% (Z-projection). 
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Figure 11. For the halo Aq-E-2, the symbols are the same as in Fig. [7] 
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Figure 12. For the halo Aq-F-2, the symbols are the same as in Fig. [7] The truncated triangle pattern in the Y-projection is due to 
naked cusps of the central caustic. The strong violation in the Z-projection is mainly caused by subhaloes with rrigub < iti^ Mq with 
a violation rate P(i?cusp > 0.187) = 6.7%. 
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Figure 13. For the halo Aq-A-2 at 2=0.6, the symbols are the same as in Fig. [7] The strong violation in the V-projection is caused 
subhaloes with ra^^^y < lO^h^^ Mq; a cusp violation rate is P(-Rcusp > 0.187) = 56%. 
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Figure 14. Effects of substructure lensing as a function of the lower cutoff subhalo mass for the halo Aq-D-2 along the Z-projection. 
The upper panels show the projected substructures with masses above a threshold. The projected centres of the subhaloes and the 
corresponding critical curves are plotted at the top. FVom the left to the right, the lower cutoff subhalo mass changes from lO*?t~^M0, 
lO'^ Mq, 10^h~^ Mq to 10^h~^ Mq. The bottom panels show the corresponding probability distribution functions of -Rcusp. Most 
substructures that survive and are projected within the central few kpc are low-mass subhaloes {< 10* h~^MQ), which dominate the 
violation of the cusp-caustic relation. 




Figure 15. Distribution of halocentric distances of the subhaloes projected within 0.05 r2oo 2.5 times the Einstein radius, indicated 
by the dotted line in each panel). The solid lines give the median spherical halocentric distances of the subhaloes; all are around 0.2 r200- 
The average number of subhaloes N is indicated inside each panel. 
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Table 3. Substructure-lensing parameters of Aquarius haloes: 



Halo Name 




13 


■S'o 


/sub,annu 


*^sub,max 


P(Rcuap > 0.187) 
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X-projection 
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5.90 


4.85 
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Y-projection 


0.657 


0.78 


0.092 


0.69 


0.059 


1.26 
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10® t 


Z-projection 
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0.76 


0.091 


0.01 


0.055 


0.08 


6.59 


10^ ^i. 


Aq-B-2 


















X-projection 


0.306 


0.73 


0.076 


0.42 


0.008 


5.91 


2.16 


— 


Y-projection 


0.408 


0.91 


0.071 


0.48 


0.007 


64.13 


0.41 
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Z-projection 


0.286 


0.64 


0.080 


0.09 


0.012 


0.09 


0.77 


— 


Aq-C-2 


















X-projection 


0.846 


0.91 


0.085 


0.13 


0.015 


19.09 


1.07 


10® ^ 


Y-projection 


0.571 


0.60 


0.107 


0.06 


0.002 


13.98 


21.41 


— 


Z-projection 


0.589 


0.70 


0.104 


0.03 


0.007 


3.71 


10.58 


10^ ij- 


Aq-D-2 


















X-projection 


0.576 


0.83 


0.101 


0.13 


0.006 


3.79 


2.13 


10^ ^ 


Y-projection 


0.655 


0.91 


0.088 


0.06 


0.005 
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0.57 


lO'' ^ 


Z-projection 
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0.79 


0.101 


0.40 


0.011 


30.58 


4.48 


10® ^ 


Aq-E-2 


















X-projection 


0.473 


0.69 


0.076 


0.18 


0.020 


3.04 


7.62 


10® ^ 


Y-projection 


0.548 


0.79 


0.056 


0.13 


0.007 


5.40 


3.29 


10^ ^ 


Z-projection 


0.474 


0.82 


0.069 


0.05 


0.006 


0.65 


1.59 


10^ ^ 


Aq-F-2 


















X-projection 


0.416 


0.86 


0.080 


0.19 


0.021 


5.83 


0.67 


10® 


Y-projection 


0.370 


0.67 


0.091 


0.09 


0.009 


1.38 
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Z-projection 


0.435 


0.85 


0.088 


0.22 


0.012 


6.74 


0.78 


10® ^ 


Aq-A-2 (Z = 0.6) 


















X-projection 


0.568 


0.71 


0.079 


0.11 


0.008 


0.60 


7.75 


10® ^ 


Y-projection 


0.731 


0.89 


0.054 


0.70 


0.022 
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Z-projection 


0.592 


0.69 


0.082 


0.33 


0.083 
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Note: Cols (2-4): bi, q-j, and S'o, the Einstein radius, axis ratio and core radius of the fitted isothermal ellipsoid (see eq. 1111 ): Col 
(5): /sub,annu is the subhalo mass fraction within a 0.1"-annulus around the outer critical curve; Col (6): asub,max is the maximum 
magnitude in the projected central 2" X 2" region of the deflection angle due to all substructures within r2oo, usually found close to an 
individual subhalo; Col (7): P(i?cusp > 0.187) is the probability (in per cent) for sources with A6 < 90° (defined as "cusp sources") 
to have iJcusp > 0.187, referred to as the "cusp-caustic violation probability"; Col (8): fa- {^9 < 90°) is the cross-section fraction (as 
defined in eq. |12) ) in the source plane for producing three close images with opening angle A9 < 90°; Col (9): A^sub.cr is the critical 
subhalo mass that causes the strong violation of the cusp-caustic relation. Arrows indicate "above" or "below" . Lenses with naked cusps 
of the caustic always have low cusp-caustic violation probability and are labelled as " — " . 
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